Gravitational Recoil from Binary Black Hole Mergers: the Close-Limit Approximation 
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The coalescence of a binary black hole system is one of the main sources of gravitational waves 
that present and future detectors will study. Apart from the energy and angular momentum that 
these waves carry, for unequal-mass binaries there is also a net flux of linear momentum that implies 
a recoil velocity of the resulting final black hole in the opposite direction. Due to the relevance of 
this phenomenon in astrophysics, in particular for galaxy merger scenarios, there have been several 
attempts to estimate the magnitude of this velocity. Since the main contribution to the recoil 
comes from the last orbit and plunge, an approximation valid at the last stage of coalescence is 
well motivated for this type of calculation. In this paper, we present a computation of the recoil 
velocity based on the close- limit approximation scheme, which gives excellent results for head-on and 
grazing collisions of black holes when compared to full numerical relativistic calculations. We obtain 
a maximum recoil velocity of ~ 64 km/s for a symmetric mass ratio n — M 1 M 2 /(M 1 + M 2 ) 2 ~ 0.19 
and an initial proper separation of 4 M, where M is the total ADM mass of the system. This 
separation is the maximum at which the close-limit approximation is expected to provide accurate 
results. Therefore, it cannot account for the contributions due to inspiral and initial merger. If we 
supplement this estimate with PN calculations up to the innermost stable circular orbit, we obtain 
a lower bound for the recoil velocity, with a maximum around 84 km/s. This is a lower bound 
because it neglects the initial merger phase. We can however obtain a rough estimate by using 
PN methods or the close-limit approximation. Since both methods are known to overestimate the 
amount of radiation, we obtain in this way an upper bound for the recoil with maxima in the range 
of 220 — 265 km/s. We also provide non- linear fits to these estimated upper and lower bounds. 
These estimates are subject to uncertainties related to issues such as the choice of initial data and 
higher effects in perturbation theory. Nonetheless, our estimates are consistent with previous results 
in the literature and suggest a narrower range of possible recoil velocities. 

PACS numbers: 04.30.Db, 95.30.Sf, 98.35.Jk, 98.62.Js 



I. INTRODUCTION 

The inspiral and merger of binary black holes sys- 
tems is one of the most interesting sources of gravita- 
tional waves that both earth-based interferometric anten- 
nas (LIGO [U, VIRGO 0, GEO600 and TAMA [|) 
and space-based ones (LISA 1) will detect. These waves 
carry both energy and momentum away from the system, 
leading to the adiabatic shrinking of the orbit, due to the 
former, and a recoil of the merged object by conserva- 
tion of momentum, due to the latter. The magnitude of 
this recoil is of astrophysical importance because it de- 
termines whether the merged hole will be ejected from 
its host galaxy. 

Possible observational evidence for such a recoil may 
be the observations of faint galaxies [y, Qj where the lack 
of a dense nucleus has been associated with the central 
black hole being ejected after merger p|. There is also 
evidence of an ejection of a supermassive black hole in on- 
going galaxy mergers, either because of recoil or because 
of slingshot due to the presence of 3 or more supermas- 
sive black holes in the merger Q . The gravitational recoil 
has also been shown to have important consequences in 
hierarchical merging scenarios and the observable struc- 
ture of galaxy nuclei. Recoil velocities of a few hundred 



km/s could be large when compared to escape velocities 
of dwarf galaxies, globular clusters and dark matter ha- 
los [1, EH • For super-massive black holes at the centers of 
galaxies, a kick of this magnitude couldpotentially trans- 
fer energy to the stars in the nucleus [8| . There are thus 
very important astrophysical aspects that can be refined 
or clarified with a better understanding of the black hole 
kicking process. 



Mass distributions without symmetries that undergo 
gravitational collapse of any sort will exhibit momen- 
tum ejection and recoil of the center of mass of the rem- 
nant due to the strong emission of gravitational radia- 
tion [H El El El El- Of particular interest is the 
case of unequal-mass black-hole binary systems. An in- 
tuitive picture of how the system gets a kick after the 
merger is the following [Trl E3| : From the center of mass 
point of view, the lighter black hole will move faster than 
the heavier one, and hence it will beam forward grav- 
itational radiation stronger. Then, there will be a net 
flux of linear momentum carried by the gravitational ra- 
diation in the direction of the lighter black hole, and 
this will cause a recoil of the center of mass in the op- 
posite direction. The first analytic studies of this sub- 
ject were carried out by Fitchett and Detweiler H 19|; 
Oohara and Nakamura 20] ; Nakamura and Haugan 21 1 ; 
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and Wiseman [T(|. Due to the strong non-linearity of 
the merger phase, analytic studies have difficulties in 
obtaining an accurate estimate of the recoil velocity. 
The first quasi-Newtonian analytic calculations were pre- 
sented in [ID, [l9j], while a post-Newtonian (PN) anal- 
ysis have been carried out in [H US, H!l- Estimates 
using black-hole perturbation theory have been given 
in [I?], [24[ , and a estimate that combines full numerical 
relativity a nd p erturbation theory, the Lazarus approach, 
is given in [251 ]. 

Full numerical relativistic simulations are a natural ap- 
proach to this problem since they can in principle handle 
the non-linearities of the gravitational field during the 
merger. The challenge is the resolution that the com- 
putational resources impose. Some calculations have al- 
ready been carried out in different scenarios to estimate 
recoil velocities. The first one was done by Anninos and 
Brandt [26| for the case of the head-on collision of two 
unequal-mass black holes. Their numerical calculations 
were effectively 2-dimensional since they made use of 
the axisymmetry of the configuration. Using the same 
type of numerical calculations they also estimated the 
gravitational radiation recoil from highly distorted black 
holes (27[ . More recently, and due to the significant ad- 
vances in 3-dimensional numerical relativity in the binary 
black hole problem [28|, [H [3(| , estimates of the radiation 
recoil velocity have also appeared [3l|, HH • 

Each of the approaches described above has its own 
limitations. Analytic approaches are able to provide ac- 
curate estimates in their region of validity. However, 
the largest contribution to the recoil velocity occurs dur- 
ing merger, precisely where the approximation methods 
break down. Numerical simulations, in principle, have 
the opportunity of producing estimates with a minimal 
number of assumptions. However, as we have mentioned, 
these calculations have also limitations and use initial 
data that is only an approximation to the actual astro- 
physical configurations. Therefore, the error bars on the 
computed distribution of recoil velocities relative to the 
distribution present in nature are believed to be large 
and, even worse, are difficult to estimate, ft is then not 
surprising to find disagreements on the estimated recoil 
velocity as calculated with different methods. 

In this paper, we present estimates of the recoil ve- 
locity using an approach not used before, the close-limit 
approximation (CLA), which combines both analytical 
and numerical techniques. The CLA was introduced by 
Price and Pullin (33[, who showed that this approxima- 
tion method provides accurate results compared to those 
obtained from numerical relativity |34l| for head-on colli- 
sions of two black holes (see also [35[ ). The CLA scheme 
is based on the assumption that in the last stage of co- 
alescence, when the black holes are sufficiently close to 
each other, the system behaves, up to a certain degree 
of approximation, as a single distorted hole. Then, the 
CLA scheme consists of establishing an appropriate cor- 
respondence between the binary black hole system and a 
single perturbed hole. Once this correspondence is made, 



one can extract initial data that can be evolved by the 
perturbative relativistic equations. From the outcome 
of the evolution, one can estimate the fluxes of energy, 
angular momentum, and linear momentum carried away 
to spatial infinity by the gravitational radiation emitted. 
The CLA scheme has been developed and applied by a 
number of authors [H, M, MM, EJ SI SI, S| . In par- 
ticular, Andrade and Price [4J| used the CLA to estimate 
the recoil velocity of a head-on collision of unequal-mass 
black holes starting from rest. 

Since the CLA applies to the last stage of the merger 
of two black holes, it is very appealing to use it to esti- 
mate the recoil velocity of the merger of an unequal-mass 
black hole binary system. With this scheme, we obtain 
a maximum recoil velocity of ~ {17,33,64} km/s for a 
symmetric mass ratio r\ = M 1 M 2 /(M 1 + M 2 ) 2 ~ 0.19 
and initial proper separations of {3, 3.5,4} M, with M 
the total ADM mass. Beyond a proper separation of 
4M the CLA is not expected to provide accurate re- 
sults [44|]. Therefore, this method cannot account for 
the contributions during the inspiral and initial merger 
phase. Supplementing this estimate with PN calculations 
up to the innermost stable circular orbit (ISCO), we ob- 
tain a lower bound for the recoil velocity, with a maxi- 
mum of ~ 84 km/s. This lower bound neglects the initial 
merger phase, for which we can obtain an approximate 
estimate by using either PN methods or the CLA. Since 
both methods are known to overestimate the amount of 
radiation during the early merger phase, we obtain, thus, 
an upper limit for the recoil with maxima in the range of 
220 — 265 km/s. We also perform non-linear fits to these 
bounds and obtain 

v nt = ai] 2 y/1 - 4ri (l + bi] + err 2 ) , (1) 

where a = 7782 km/s, b = -2.507 and c = 2.727 for 
the lower bound and a = 14802 km/s, b = —1.1339 and 
c = 1.4766 for the upper bound. 

The plan of this paper is as follows: Sec. |TT] describes 
the main procedure of our calculation; Sec. [Hi] constructs 
initial data for a quasicircular binary black hole system in 
the 3 + 1-formalism; Sec. IIVI maps these initial data to a 
single perturbed black hole spacetime, such that it is suit- 
able for a CLA evolution; Sec. IVl describes the numerical 
implementation and presents results from the evolution 
within the CLA scheme; Sec. I VII estimates the lower and 
upper bounds, as well as constructing the non- linear fits 
to these bounds; we finish in Sec. IVIII with a summary 
and discussion of the main results and points to future 
research to obtain improved estimates. 

The conventions that we use throughout this work are 
the following: For the 4-dimensional spacetime, we use 
Greek letters for the indices and a semicolon for the co- 
variant derivative. The Schwarzschild metric admits a 
2+2 decomposition consisting of the warped product of 
a Lorentzian 2-dimcnsional manifold with the 2-sphere 
(see 0, S|). O n the 2-dimensional Lorentzian manifold 
indices are denoted by capital Latin letters, the covari- 
ant derivative associated with the 2-dimensional metric 
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is represented by a vertical bar, and the Levi-Civita an- 
tisymmetric tensor by e AB . On the 2-sphere indices are 
denoted by the lower-case Latin letters a,b,...,h, the co- 
variant derivate by a colon, and the Levi-Civita antisym- 
metric tensor by e ab . When using the 3+1 decomposition 
of spacctime quantities, spatial indices are denoted by the 
lower-case Latin letters i,j,k,.... Uncontrolled remain- 
ders are denoted with O(A) or 0(A, B), which stands for 
terms of order A and terms of order A or B respectively. 
Although usually, when dealing with order symbols, A 
and B must be dimensionless, here they will not be, but 
can be made to be dimensionless through division by the 
appropriate factor. We also use physical units in which 
G = c = 1. 



II. DESCRIPTION OF OUR CALCULATION 

In this paper, we use the CLA scheme to calculate the 
recoil velocities after the merger of an unequal-mass bi- 
nary black hole system. Due to the complexity of the 
calculation, we discuss here the different steps involved, 
while getting a glimpse of the general scheme. First, 
we need to construct initial data corresponding to a non- 
spinning binary black hole system in quasicircular orbital 
motion. The method employed to construct the data is 
the standard one: we solve the constraints on an initial 
slice using the York-Lichnerowicz conformal decomposi- 
tion. Then, the solution needs to be expanded in two 
parameters: the separation of the two black holes, based 
on the main assumption of the CLA, that is, small separa- 
tion between the holes; and their linear momenta, rooted 
in an additional slow motion approximation (47j . 

The second step is to establish a map between this 
initial data and the generic initial data corresponding 
to a perturbed single black hole. In this work we only 
consider the case in which the single black hole is a non- 
rotating Schwarzschild hole. There is also the possibility 
of considering a Kerr black hole (see [53] for details), 
but the CLA machinery in that case is more intricate. 
After expanding the initial data in the separation and 
linear momenta, it is straightforward, after some coor- 
dinate changes, to identify a Schwarzschild background 
and its perturbations. 

Once the perturbations have been identified, we need 
to calculate initial data suitable for evolving the lin- 
earized (around the Schwarzschild background) Einstein 
equations. The spherical symmetry of the background 
allows us to separate the linearized equations. Then, by 
decomposing the perturbations in spherical harmonics we 
obtain decoupled equations for each mode. Moreover, by 
appropriately reparameterizing the perturbations, we can 
decouple the equations for each individual mode, so that 
the problem reduces to solving a master equation for a 
complex combination of the metric perturbations. These 
master equations (usually known as the Regge- Wheeler 
and Zerilli-Moncrief equations) are 1-dimensional wave- 
type equations containing a potential that accounts for 



the effect of the background spacetime curvature. There- 
fore, the problem of providing initial data reduces to find- 
ing initial conditions for these master functions. 

The initial data contains three parameters that we 
need to specify. These parameters are associated with the 
initial distance between the holes, the mass ratio, and the 
initial linear momentum. The mass ratio is an indepen- 
dent parameter that will be used to study the functional 
behavior of the recoil velocity. The distance and linear 
momentum parameters determine the dynamical charac- 
ter of the binary and, therefore, they must be chosen 
carefully. To that end, we use the standard method of 
minimizing the binding energy of the system, so that the 
binary is in a quasi-circular orbit. The expressions that 
we obtain are formally the same as in Newtonian theory, 
although they cannot be assigned the same interpreta- 
tion, since they are expressed in terms of bare parame- 
ters. In order to relate these parameters to meaningful 
physical ones, we must introduce a proper separation and 
a physical mass ratio. The proper separation can be cal- 
culated by evaluating the minimum proper distance be- 
tween the marginally trapped surfaces surrounding each 
individual hole. 

The final step is to solve the master equations and 
evaluate the different physical quantities of interest. The 
metric waveforms h + and h x , together with the fluxes 
of energy, angular and linear momentum carried away 
by gravitational waves can be computed in terms of the 
master functions and their first time derivatives. In this 
paper, we include the general formulae for the linear 
momentum fluxes in terms of the perturbation master 
functions. We present several plots of these quantities, 
together with plots of the recoil velocities for different 
initial separations. 

III. INITIAL DATA 

In this section, we begin the initial data construction 
for an unequal-mass binary black hole system suitable 
to the CLA scheme. To that end, we extend the results 
of Andrade and Price (44[, who carried out the calcu- 
lation for unboosted head-on collisions, and also extend 
the results of Khanna et al [47|, who constructed data 
for equal-mass black holes in a quasicircular orbit. Our 
calculation not only allows for arbitrary mass ratios, but 
it also includes higher-order terms in the expansion of 
the initial data, which are essential in the calculation of 
the recoil. 

In order to solve the Hamiltonian and momentum 
constraints, we use the conformal transverse-traceless 
method of Lichnerowicz, York and others [48l. |4^. IHol. l5ll. 
l52j . The 3-metric 7^ is decomposed in terms of a con- 
formal factor <!> and an auxiliary metric % • , 7^ = 3> 4 7jj , 
which here we assume to be conformally fiat: 

ds 2 = j ij dx i dx j = $ 4 (di? 2 + R 2 dil 2 ) , (2) 
where dfl 2 = dO 2 + sin 2 6d<p 2 is the line element of the 
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2-sphere. For the extrinsic curvature , we choose a 
maximal initial slice, that is, is trace free: j^K^ — 
0. Then, we also conformally decompose the trace-free 
extrinsic curvature K 4 a as 



(3) 



and we further make the choice that the longitudinal part 
of K { j vanishes, so that is a symmetric transverse 
traceless tensor. Then, the momentum and Hamiltonian 
constraints reduce to 



0. 



(4) 
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where V i and V 2 denote the covariant derivative and 
Laplacian associated with the flat 3-metric ■ . The mo- 
mentum constraint [Eq. Q] can be exactly solved using 
the method of Bowen and York 53]. For a single black 
hole located at R = R Q with linear momentum P it can 
be written as follows: 



j>one _ 

« ~2|i2-ilJ2 



(6) 

P' is the ADM momentum of a single hole, while n z 
is a unit vector in flat three-dimensional space directed 
from the location of the single hole to an arbitrary point, 
namely 



\R-R n 



(7) 



and the vertical bars, | • | , denote the norm of vector in the 
flat 3-dimensional space. In order to construct a solution 
for two holes, we can simply superpose two solutions of 
the type of Eq. ©. 

Before constructing the extrinsic curvature, it will be 
uscfull to first describe the initial physical configuration. 
The system we are modeling consists of two black holes 
with masses M x and M 2 located on the X-axis, a coor- 
dinate distance d apart, as shown in Figure [TJ In this 
figure, R 1 , R 2 , and R are radial vectors that point from 
the origin to hole 1, hole 2, and an arbitrary point, re- 
spectively. Moreover, R 12 = R 2 — R\ is also a vector 
that points from hole 1 to 2, and P and — P are the lin- 
ear momenta associated with holes 1 and 2 respectively. 
Since the linear momenta are parallel to the Y-axis, the 
orbital angular momentum is directed along the Z-axis. 
With such physical scenario, the solution of Eq. (H|) can 
be written as (see also (47j): 



K, 



K^[R = R^P] 



R 2 ,-P], (8) 



where we have defined P = Py. The ADM momentum 
corresponding to K°™ e is simply P and the one associated 

with K11 is zero. 



FIG. 1: Schematic diagram of the initial physical configura- 
tion. The linear momenta are parallel to the Y-axis and span 
the X-Y plane, so that the angular momentum is aligned with 
the Z-axis. 



Let us now concentrate on the solution to the Hamil- 
tonian constraint. Using Eq. |8]) in the Hamiltonian con- 
straint leads to a source term quadratic in P. Wc now in- 
troduce a slow-motion approximation, where we assume 
that the linear momentum P is small, in the sense that 
terms of 0(v 2 ) are much smaller than terms of 0(v), 
where v is a measure of the orbital velocity. We, thus, 
neglect terms of 0(P 2 ), so that the Hamiltonian con- 
straint reduces to the Laplace equation 



= 0. 



(9) 



The solution of this equation depends on our choice of 
topology. If we choose the initial slice to have a single 
asymptotically flat region, the solution to the conformal 
factor is the Misner solution 54], but if one chooses the 
slice to have three asymptotically-flat regions, the solu- 
tion is the Brill-Lindquist (55j one. In this paper, we 
adopt the latter and the conformal factor takes the form 
of the Newtonian potential: 



$ = 1 



2\R-R 1 \ 2\R R 2 



(10) 



where m 1 and m 2 denote the bare masses of each indi- 
vidual hole. One reason for choosing the Brill-Lindquist 
(BL) solution is that it is simpler to manipulate, while it 
has also been shown [44[ to lead to essentially the same 
results when calculating recoil velocities for head-on col- 
lisions. We remark that, although terms of C(P 2 ) have 
been neglected, they can be straightforwardly added in a 
perturbative fashion, but this will be studied elsewhere. 

Let us comment further on the topology of the initial 
slice associated with the BL solution, as it is important 
in some calculations. As we have already mentioned, 
this solution has three asymptotically-flat regions: one 
of them, S , corresponds to the region far from the two 
holes, R = \R\ > \R X \ = R 1 and R > |H 2 | = P 2 ; the 
other two, H 1 and E 2 , are associated with hole 1 and 
2 respectively. By simple inspection of the conformal 
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factor [Eq. (fT0|) ]. the solution seems ill-behaved at the 
location of the holes, R — R x and R = R 2 , although it 
is actually not. Near each hole, the geometry is invari- 
ant under the transformation: \R— R&\ — > m A /(4R' A ) 
(A = 1,2). The value m A /2 coincides with the inter- 
section of the event horizon with the initial slice for a 
single hole and it is a fixed point in the transforma- 
tion. This value is sometimes referred to as the throat, 
joining two asymptotically-flat regions. Therefore, the 
points R = R 1 and R = R 2 are just an image of the 
infinities of T, 1 and E 2 . For a single hole, there are two 
asymptotically- flat regions, and its mass, equal to m, is 
the same independent of which region we evaluate it on. 
In the case of a binary system, the gravitational inter- 
action between the holes will change the value of the 
individual masses. Actually, there is not an invariant 
measure of them in S , but such a measure does exist 
on Ej and £ 2 . Doing the calculation yields the following 
result [55| 

M, = m x (l + ^) , M 2 = m 2 (l + ^) , (11) 



2d 



2d 



where d = |-Ri 2 | • In £ i we can compute the total mass 
of the binary system, the ADM mass of the system. We 
call the result M and it is given by 



M = m 1 



(12) 



Eqs. ©, (Unj), and © are the initial data. We should 
note that, apart from our choice of initial data (BL 
conformal factor and Bowen-York extrinsic curvature), 
there are other possible choices that can be used in the 
CLA scheme. Some examples of other possible data 
sets are the following: a Misner conformal factor with 
a Bowen-York extrinsic curvature with inversion symme- 
try through the throats; Kerr-Schild initial data [42 ] . 

The next step in the construction of the initial data is 
to put it in a form suitable for the CLA scheme. Before 



doing so, however, it is convenient to study the param- 
eters that determine the configuration described by the 
data. To begin with, let us introduce the bare mass ratio 



(13) 



The initial configuration can then be fully specified in 
terms of the parameters (g, d, P) . Since q and d are 
bare parameters, in the sense that we cannot give them 
a physical meaning, we are going to introduce analogous 
parameters, which can be given a physical interpretation. 
First, we introduce the mass ratio between the individ- 
ual masses of the holes as computed in their respective 
asymptotically-flat regions: 

The quantities q and Q are related through d via [44j ]: 
l+[(l + g)/2](M/d) 



Q = 



Defining a useful distance is a more difficult matter. In 
this paper, we use the same distance as Andrade and 
Price [44j, which is the proper distance between the 
points where the marginally trapped surfaces, surround- 
ing each individual hole, cross the X-axis (we obviously 
refer to crossing points closer to the opposite hole). When 
the initial configuration does not present a common ap- 
parent horizon, these marginally trapped surfaces are the 
individual components of the apparent horizon. If x 1 and 
x 2 stand for these crossing points, the distance we have 
just defined is given by D — J x 1 $ 2 dx , which yields 




in 1 m 2 
4(X 1 -x 1 )(X 1 -x 2 ) + 4( Xl -X 2 )(x 2 



where Ax = x x —x 2 > , and X x (= i? 1 ) and X 2 (= — R 2 ) 
are the x-coordinate locations of holes 1 and 2 respec- 
tively in the conformal space. In summary, we can de- 
termine our initial configuration either by specifying the 
set (q, d, P) or the set (Q, D, P) . 

The CLA scheme assumes that the black holes are suf- 
ficiently close enough, which allows us to expand the ini- 
tial data in R S> R x and R 3> R 2 . In this sense, it is 
useful to choose the coordinate origin, R = , in such a 
way that it coincides with the bare center of mass of the 



1 

binary black hole, that is: 

m 1 R 1 +m 2 R 2 = 0. (18) 

We can write then R x and R 2 in terms of the separation 
vector R 12 via 

R 1 =XiR 12 , -R 2 = x 2 -R 12 , (19) 
where we have defined 
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Then, it is natural to expand the initial data in the fol- 
lowing dimensionless parameter e 



e= € = 



R 



(21) 



A key formula for performing these expansions is the fol- 
lowing: 



(22) 



£=0 



where we have introduced the following unit vectors 

R = f , e=-, (23) 



and where denote the Gegenbauer polynomials. 

These polynomials, also known as ultra-spherical polyno- 
mials, are a generalization of the Legendre polynomials 
to {N/2 + 2)-dimensional spaces, which are common in 
angular momentum theory [5a . |57| . For the special cases 
of N = 0,1,2, these polynomials reduce to Legendre 
polynomials and Chebyshev polynomials of type 1 and 2 
respectively. We refer to Appendix lAl for more details on 
these polynomials. 



We now use all these definitions and results to expand 
the conformal extrinsic curvature given by Eq. (jHJ) in e to 
arbitrary order. The result we obtain can be written as 
follows: 



q °° / J \ l + 1 

k = 2^e(1) (/i +i "xn{c[2w-( p ^)i,]+^( p ^)M 



Cf ,2) [{P-e)%-2P {i e j) 



-C. 



(5/2) 



(P-e) R i R j +2{P-R)R {i e j) ]} 



1+2 



1+3 



X{ +2 ~ Xi +2 ) Cf' 2) {P ■ R) Mj + 2{P ■ i) R 



(i e i) 



{x{ + "-xr)cf /2 \P-e)e l i 1 , 



(24) 



where for simplicity we have omitted the argument of the 
Gegenbauer polynomials, which still is e • R. It is worth 
noting that the e° term has identically vanished due to 
the choice of coordinate origin, which coincides with the 
bare center of mass. Another interesting fact is that only 
combinations of the (3/2)- and (5/2)-Gegenbauer poly- 
nomials appear due to the combination of odd powers in 
the denominators of the extrinsic curvature. 



In this paper we are going to consider terms up to 
0(e 3 ), which is enough to get a gravitational recoil effect 
and actually the dominant part of it (see, e.g. [3). Ex- 
tensions of our calculations to higher order are in prin- 
ciple straightforward, but we are not going to present 
them here. Then, up to this order of approximation we 
can rewrite Eq. (|24|) as follows: 



K, 



9 d 
2^3 
3 d 
2^3 
15 d 

y^3 

3 d 
2R3 



Re + lH-^-[5{R-e) 2 -l]} (P •/?!-,, 2P (l R 3) 



1 + 3 



Re 

1 + 5 



2i?l + q 

^—iR-i 

Rl + q 

1 d 1 -q 



2Rl + q 
dl-q A 



(P • e)% - 2P {i e j} 
[7{R-i) 2 -l}){P-R)RA 



Rl + q 



-Re 



3 d 2 1 - q 
2R 4 l + q 



(P • e)R,R 3 + 2(P • R)R (i e j) 

(P • R)e t e 3 + 2(P • e)^] + 0{Pd 3 ) . 

I 



(25) 



The lowest-order contribution is of 0{Pd) and it is the only contribution used by Khanna et al [47f for grazing 
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collisions of equal-mass black holes. The next contribu- 
tion is of 0(Pd 2 ) and, as far as we know, this is the first 
time it has been considered. 

Let us now look at the conformal factor [Eq. I|10p]. 
Using Eq. (|22|) . we can also expand $ in Gcgenbauer 
polynomials to obtain 

* = 1 + H + E ^ 1/2) ' g ) e ' M + -2X0 -(26) 



£>2 



It is important to recall that in solving the Hamiltonian 
constraint we have used a slow-motion approximation, 
neglecting terms of 0(P 2 ). The terms we are, thus, ne- 
glecting are of 0(P 2 d 2 ). This expression can also be 
written in terms of the parameters (q, d) and M, and in 
terms of Legendre polynomials. In this way we obtain 



$ = 1- 



M 
2R 



e>2 



P t (R-e) + 0{P z d 2 ), (27) 



where Pe denotes the Legendre polynomials and where 
the coefficients <p e are given by 



1 



{(-lY + q"- 1 } 



(1 + 9)' 



(28) 



The I = 1-term vanishes due to the choice of the origin 
of coordinates in Eq. (fT8|) . Finally, the expansion of the 
conformal factor up to third order in d is given by 



1 

$ = 1 + — + - 



M 
2R 



q 



Md 7 



\q{\ 



2(1 + q) 2 
q) Md 3 



2(l + q) 3 i? 4 



R 3 

Ps(R-e) 



0(P 2 d\d 4 



(29) 



With this, we finish the construction of initial data to be 
used in the CLA scheme. To summarize, we remark that 
this construction is based on expansions on two different 
parameters: P (related to the slow-motion approxima- 
tion) and d (related to the assumption that the holes are 
close to each other). Since P and d have dimensions, the 
meaning of these expansions is that terms of order d N 
and/or P M are smaller than terms of order d N_1 and/or 
P^ 1 ^ 1 . As we are going to see later, these expansions will 
provide the leading contribution of the multipoles £ — 2 
and £ = 3. 



IV. THE CLOSE LIMIT APPROXIMATION 

The next stage in our computation is to recast the 
initial data just constructed into data for a perturbed 
Schwarzschild black hole, which is the essence of the 
CLA scheme. In this way we can extract initial data to 
be evolved by the corresponding perturbation equations. 
Thanks to the expansions performed in the previous sec- 
tion, the main task now becomes the extraction of the 
different multipoles from the data. 



The 3-metric on the initial slice is conformally flat and 
hence determined by the conformal factor <I>. If we look 
at the lowest-order contribution [see Eq. (|2"9"]) ] we real- 
ize that it coincides with the 3-metric of Schwarzschild 
spacetime associated with the {t — const. }-slicing in 
isotropic coordinates, being t the Schwarzschild time co- 
ordinate. However, in order to make the connection with 
perturbation theory, it is very convenient to reexpress the 
initial data in Schwarzschild coordinates: 



ds 2 = f^dr 2 + rW , / 



2M , x 

1 , 30 

r 



where we recall that M is the total ADM mass. 
The transformation from isotropic coordinates to 
Schwarzschild coordinates is given by the following re- 
lations: 



1 



R 1 



1 1 • <*> 



Applying this transformation to the 3-metric of our initial 
data we obtain: 



ds 2 



T A {r l dr 2 + rW) 



where 



T = 



1 + M. 

1 T 2R 



(32) 



(33) 



In order to construct initial data for the perturbations, 
evolve it, and compute from the result all the relevant 
physical information, in the next subsections we give a 
summary of (non-rotating) black-hole perturbation the- 
ory and the main tools needed for the application of the 
CLA scheme. Afterwards, we apply this machinery to 
the construction of the initial data and describe how the 
energy, angular momentum, and linear momentum fluxes 
carried away by the gravitational waves are evaluated. 



A. Black hole perturbation theory 

The CLA is based on the fact that, in the last stages of 
coalescence, the gravitational field can be modeled, to a 
good degree of approximation, as the gravitational field 
of a single perturbed black hole. Thus, perturbation the- 
ory plays a key role in our calculations and it is worth 
reviewing its main concepts and tools. The starting point 
is the assumption that the spacetime metric, g^ u , can be 
written as: g = g^ 1 + h v , where g^J denotes the 
background Schwarzschild metric and the first-order 
perturbations. Then, we can take advantage of the spher- 
ical symmetry of the Schwarzschild metric to simplify the 
structure of the perturbations and of the equations that 
govern them. We can do this by expanding the pertur- 
bations in tensor spherical harmonics. It turns out that 
the linearized Einstein equations (in this case, around the 
Schwarzschild background) decouple for each harmonic. 
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Not only this, we can distinguish between the pertur- 
bative modes with polar parity, which pick up a factor 
of (—1)' under parity transformations, and the ones that 
have axial parity, which pick up a factor of (— This 
distinction is important because polar and axial modes 
also decouple. 

Following this discussion, we split the metric pertur- 
bations h^y into polar and axial perturbations, h^ v = 
/i* + h^ v . And these perturbations can be expanded in 
tensor spherical harmonics as 



where 



l.m 



£m i p _ i p,£m 



hf 1 S e a m 



tra otra 



(34) 



(35) 



where v A m 



£m 

Pa 



l^Cm 



(r 2 /2)G tm 



1 Trim 

- 2 H\A 



And for axial modes 



(42) 



The equations for the metric perturbations decouple in 
terms of complex master functions, so that once we solve 
the decoupled equations for these master functions all the 
metric perturbations can be reconstructed from them. 
In the case of axial modes, it was first done by Regge 
and Wheeler [58[ , and for polar modes by Zerilli [5f| and 
later by Moncrief [6fJ. These functions are made out of 
metric perturbations and their first derivatives and they 
are gauge invariant. It is also possible to express them in 
a covariant form. In the case of polar modes, the Zcrilli- 
Moncrief function can be written as follows \6M 



Ira 



ab 



1 + A, 



K 



I III 



1 

T. 



r n AB 



^\A 



(,43) 



ip,<m 

''(HI/ 



u£m \r£m 
n AB Y 



Pa Y a 



1 ab 



,(36) 



where asterisks are used to denote components that are 
given by the symmetry of these tensors. Y em are the 
scalar spherical harmonics [see Appendix I A 21 for the con- 
ventions that we use and other details] . Y^ and S^ 71 are 
vector spherical harmonics and are defined (for I > 1) in 
terms of the scalar spherical harmonics by 



Y. 



Im 



Y. 



s, 



Ira 



b-w-tm 
i *b 



(37) 



Finally, Y* 



and Sf™ are (symmetric) tensor 



spherical harmonics, which can also be defined {Z^™ and 
gim on iy f or i > 2) in terms of the scalar spherical har- 
monics by 



\f£m 
1 ab 



Y em n n 



z. 



ab 



\rira 
1 :ab 



£(£+!) 



Y" m n ab ,(3S) 



where X e = (£+2)(£-l)/2 and A e = \ e + 3M/r. For axial 
modes, instead of using the well-known Regge- Wheeler 
master function 



= -L r \ A h 



\A : 



(44) 



we are going to use the master function introduced by 
Cunningham, Price and Moncrief [621 ], in the form used 
in [6l|, [63| . The main reason for this choice is that it is 
simpler to evaluate the fluxes of energy, angular momen- 
tum, and linear momentum. Moreover, the contributions 
of axial modes to these physical quantities have the same 
form as the one of polar modes. Nevertheless, for the 
sake of completeness, we provide formulae for both mas- 
ter functions. The Cunningham-Price-Moncrief master 
function can be written in covariant form as 16 ill 



Cm 



AB 



A, 



u£m 
n B\A 



\A n B 



(45) 



5im 
nh 



Qirn 



(39) 



Here, the sign convention for the Levi-Civita tensor as- 
sociated with the metric of the 2-sphere is: e 6ip = sin#. 
In Appendix I A 21 we give the orthogonality relations for 
the different harmonic objects. All perturbative quanti- 
ties, scalar (hjfij), vectorial (p^ an d <Ta )? and tensorial 
(K lm , G lm , and g| m ), are functions of t and r only. 

The metric perturbations are in general not invariant 
under transformations of the mapping between the back- 
ground and perturbed spacetimes, or in other words, they 
are in general not invariant under gauge transformations. 
However, for the case of a spherically-symmetric back- 
ground, like the Schwarzschild metric, there is a complete 
set of perturbative quantities that are gauge invariant. 
For polar modes this set can be chosen as follows 



lira 
n AB 



t£m 
n AB 



ZV A\B ■ 



K 



I in 



K 



£ra + 1) q£ 



r \A 

2—vT 



(40) 
(41) 



In Schwarzschild coordinates these functions take the fol- 
lowing form (the connection with the Regge- Wheeler pa- 
rameterization of the perturbations is given Appendix [B| 



Cm 



1 + A 

I 
A, 



fh 



{K im + (l + X e )G im 
2 



< III 



rd, r K 



km 



(l + X e ) P : 



I'm 



^£m _f_ ( rim _ 1 Q rr£ 

*rw r \ 2 



r 



A, 



K m - d r h t 



Cm 



£m 



(46) 



(47) 



(48) 



These master functions obey the following wave-type 
equation with a potential: 



lm , =0 (49) 
cpm/zm u ' v^ j ; 
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where r t is the so-called tortoise coordinate (r, = r + 
2Mln(r/(2M) - 1)). The potential for the axial modes 
is the Regge- Wheeler potential 



'(r) = £ [1(1 + 1) 



r / 



(50) 



and the one for polar modes is the Zcrilli potential 



TrZM 



(r) 



/ 



2 A 2 



3M 



M 2 



71/ 



2AJ ( 1 + A* + ] + 18— I A £ + — 

r 

Once the different master functions have been com- 
puted we can estimate the energy and angular momen- 
tum carried out by the radiation field to infinity. We 
can do this by using the expressions of the energy and 
angular momentum fluxes at infinity obtained from the 
Isaacson's averaged energy-momentum tensor for gravi- 
tational waves [64|,[65[ (see also [H, H?}). In terms of the 
axial and polar master functions the expressions are 



E = 
Gw 64tt 



- E 



■2)1 



i>2,m 



(£-2)\ 



im |2 

CPM I 



I* 



(52) 



- E 



64tt 



l>2.: 



(1 + 2) 
\i-2) 



1 ( m tm W im 4- f' m f (m ) (53) 

I V CPM * CPM T ^ZM 1 ZM J V""/ 



We can also construct the metric waveforms by using 



h, — ih v 
+ x 2r 



e>2,m 



(e-2)\ 



(54) 



where _ 2 Y lm denotes the spherical harmonics of spin 



weight —2 (see, e.g. [68[ and Appendix IA 31 for details). 
In this work we are interested in studying the gravita- 
tional recoil due to the merger of unequal-mass black- hole 
binary systems and therefore, we want to evaluate the 
flux of linear momentum emitted in gravitational waves. 
This quantity can also be computed from the Isaacson's 
energy-momentum tensor and can be written in terms of 
the metric waveforms as follows: 



pk 



16tt 



(55) 



where r k obs is a unit vector that points from the source 
to the observer. We can then express the components of 



o/>-s 



in terms of scalar spherical harmonics as 



obs 



= -2W-lMR(Y^),3(YM),- 



Y 1 



V2 



(56) 



where 3? and 9 denote the real and imaginary parts of 
a complex number. By simple inspection of the linear 
momentum flux in Eq. (|55p . and taking into account the 
harmonic structure of the metric waveforms in Eq. (154 



and of f^ bs in Eq. ([56")) , we realize that all terms in the flux 
contain the product of three spherical harmonic objects. 
Therefore, in order to obtain a practical expression for 
P' w we need to use the machinery for studying coupled 
angular momenta common in quantum physics 1561 . l57j | . 
The calculation goes along the lines described in [63] , and 
some details are given in Appendix [Cj The result can be 
written in the following form 



" ^ E V + 2 W~ !) { V(* - m)(t + m + 1) (*t%M +1 - ^m^z'm +1 ) 

t>2,m 

+ ^(£ + m)(e-m + l) (^™* £ c ' P m M - 1 - *cpm*z'm - 1 ) } , (57) 



GW MlT ekm V 2 ) ! V + lV(2* + 3)(2£+l) I V 1 )[ J ^ ZM ZM CPH CPM ) 

+ y /(£- m + 2 )(e-m+l) (viZHt 1 '™- 1 + *S? M ^i m_1 ) } 
" ei^ E ( £ + 2 )( £ - !) { V 7 ^ - m)(e + m + 1) - * C pm*zm +1 ) 

f>2,m 

- y/lt + m)(e-m + l) (^MZ- 1 ~ *cPM*' M ra - 1 ) } , (58) 
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p z = 



J_ v (1 + 3)! (e + m + W-m+1) f, im ^ +1 , m 
32tt ^ (£- 2)1 y (2^ + 3)(2£ + !)(£ + I) 2 V ZM ZM 



^CPM ^ CPM 



(59) 



l>1,m 



In conclusion, all we need to extract relevant physical 
information is the master functions. In the next subsec- 
tions, we extract initial data for these master functions. 

B. Relation between ADM variables and metric 
perturbation 

In section Hill we constructed initial data for a binary 
black hole system in coalescence. The procedure used 
for this construction was based on the 3 + 1 ADM for- 
malism [69| and, hence, the initial data is given in terms 
of ADM variables [Eqs. (2$, (35J and ©]. Then, in 
order to build initial data for the evolution of the mas- 
ter functions, we need to first find the relation between 
the ADM variables and the metric perturbations (see, 
e.g. |z3|). This means that we need to use the rela- 
tions between the components of the 3-metric 7^ and the 
metric perturbations, and also the relations between the 
metric perturbations and their first derivatives and the 
components of the extrinsic curvature 1£+ . For the for- 
mer, we use the fact that the components of the 3-metric 
are the spatial components of the orthogonal projection 
operator on the hypersurfaces of the spacetime slicing, 
described by a normal n^: 

7 M „ = g» u + ^n v ■ (60) 

Then, the different modes of the harmonically decom- 
posed 3-metric are related to the metric perturbations 



via 





= -/a - , 


(61) 


it? 


— 'hr 1 ' 


(62) 




flm\rim , uCm oCrn 
— Pt X a + n t D a ' 


(63) 


V™ 
Irr 


= r l + h l ™Y lm , 


(64) 


ira 


ILmvtra , ulrn ultra 
— Pr Y a n r J 


(65) 






t~ihn ry£m\ 
r <-f ^ab ) 




+ H b ab . 


(66) 



The first three equations are related to the choice of slic- 
ing, that is, to the choice of shift vector (3 l and lapse a . 
Actually, the lapse and shift at first order are given by 

a 2 = f(l-fh tt ) = f-h tt + 0(h 2 ), (67) 
13'' = h t l + 0{h 2 ). (68) 

The relation between the extrinsic curvature and the 
metric perturbations can be found through the relation 
between the 3-metric and the extrinsic curvature 

^liV — 2°^n7piy i (69) 

where the symbol £ denotes Lie differentiation, and 
Eq. (|6H)) between the 3-metric and the spacetime met- 
ric. In this way, we find that the different modes of the 
harmonically decomposed extrinsic curvature are related 
to the metric perturbations by the following expressions: 



Cm 



K 



jy-Cm 
ra 



K ab 



\ira 



2V7 L 

r 2 f 

V71 



2h 



Iral 



f 

r 

2 



Ira 



p e r - p £ r' + -pi 

r 



Cm 



£(£ + 1) 



-Pt 



Ira 



~ylm 

l^lra 
'Hr 

2/, 



Y. 



Cm 



hi m -hi m ' + -„, 



i in 



s, 



i III 



1 ab 



'-lira 



G 



-v tm 



rytra 



2h 



Cm 



5£m 
ab 



(70) 
(71) 
(72) 



where the dots and primes denote partial differentiation C. Initial data for the metric perturbations 

with respect to time t and radial coordinate r respec- 
tively. 

Before computing initial data for the master functions, 
we must find data for the metric perturbations, that is, 
find {hny, h„ v ), in the parameterization given in Eqs. (|35j) 
and ([36)1 on the initial slice t = t a . To begin with, since 
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our 3-metric is conformally flat, the following metric per- 
turbations vanish on the initial slide: 



We have also seen that the conformal factor, the physical 
and the conformal extrinsic curvatures can be formally 
expanded in powers of d and P as follows 

* = *(o) + ^{2)d 2 + $ ( 3)rf 3 + 0{P 2 d\ rf 4 ) , (74) 
= PdK {l)l] +Pd 2 K {2)l] +0(Pd 3 ), (75) 

where and K, 2 \ij are the coefficients of the terms 

of order Pd and Pd 2 respectively in the expansion of 
K t j . Then, the physical extrinsic curvature, Ky, given 
by Eq. ([3]), can be formally expanded in the form 

K l3 = <P- 2 (Pd K {m + Pd 2 K (2)ij ) + 0(Pd 3 ) , (76) 

which means that in order to obtain the physical extrinsic 
curvature up to 0(Pd 2 ) we only need the zero-th order 
piece of the conformal factor. The explicit expressions of 
the coefficients of these expansions are given by equations 
(22), (HSJ and ©. 

With this in mind, we are going to extract the remain- 
ing modes of the initial data. To that end, we use the 
expression of the separation vector e in spherical coordi- 



nates, namely 



sm tt cos ip , ■ 



cos 6 cos ip 
~R 



sm ip 
Rcosd 



(77) 



Then, the non-zero components of the 3-metric on the 
initial slice are given by 



q Md 2 1 
1 + 2 - - — -P 2 (0 



(1 + q) 2 r 3 a 5 
q(l-q)Md 3 1 



[1 + q) 3 r 4 a 



+ 0(P 2 cf,cf), (78) 



-fab = r 2 Q ab 



q Md 2 1 
1 + 2 - q Nn — -P 2 (H) 



T 1 



[1 + q) 2 r 3 a 5 ' 

+ 0{P 2 d 2 ,d i ), (79) 



(1 + q) 3 r 4 a 
where we have introduced the following definitions 

4 = sm cos (p, a = . 



(80) 



We can now rewrite the 3-metric in terms of spherical 
harmonics as follows 



= / 



1-2*/? 



Md 2 1 



5 (1 + q) 2 r 3 
V5 3?(r 3 < 3 )]} +0{P 2 d 2 ,d 4 ) 



Y 2 ' - V6K(Y 2 ' 2 ) 



-2W?^ 



q)Md 3 1 



7 (1 + q) 3 



Vs^Iy 3 - 1 ) 



(81) 



lab = r 2 n ab -2J- 



Md 2 1 



5 (1 + q) 2 r a 5 
- ^( Y !f)} +0{P 2 d 2 ,d i ). 



Y 



2.0 



VeniY 2 , 2 ) 



_, } r ¥q(l-q)Md 3 1 



7 (1 + q) 3 r 2 



V3n(Y a f) 



(82) 



In order to repeat this procedure with the extrinsic cur- 
vature, we first need to compute the components of the 
conformal extrinsic curvature with the separation vector 



e of Eq. ([77)) . The components of the conformal extrinsic 
curvature are given by 
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Pd l — o Pd 2 

K RR = 3-^3 sin 2 gsin(2y))-3 1 + ^ Ri sin y sin (5 sin 2 6>cos 2 ip - 2) + 0(Pd 3 ) , (83) 

3 1 — a Pd 2 

K m = -— | cos 6» sin 9? (5 sin 2 6» cos 2 93 + 3) + 0(Pd 3 ) , (84) 
Pd 3 1 — Pd 2 

K Rv = sin 2 6> + - 1 - q R3 sin 6 cosy [sin 2 (5 cos 2 tp - 14) +3] + 0(Pd 3 ) , (85) 

3 Pd 3 1 — o Pd 2 

K ee = -— sm(2ip)[cos{29)- 5}- -j^^smdsmip [5 cos 2 if (cos 2 d- 3) +3]+ 0{Pd 3 ), (86) 

3 Pd 3 1 — Pd 2 
K 6(p = --— sin(2(9)cos(2y) + -- -— - sin 2 (9 cos cos <p (5 cos 2 tp - 2) + ©(Pd 3 ) , (87) 

3 Pd 1 5 1 — o Pd 2 

K vv = -— sin 2 6 sin(2y>) [1 + 3 cos(26>)] + — fl2 sin 3 flsiny [cos 2 y (3 sin 2 - 2) - l] + 0(Pd 3 ). (88) 



I 

Here we have checked that the terms of O(Pd) agree with extrinsic curvature in terms of spherical harmonies and 
those in [47j (up to a typo in their value of the {9, tp} com- Schwarzschild coordinates. This quantity is given by 
ponent). The next step is the calculation of the physical 



K„ 



V 5 r" 3 / 
2V3^4^° 1 



3tt 1 - q Pd 2 1 



7 1 + 9 r A fa 2 
1 FTl-qPd 2 1 



VT43(y M ) + ^{Y 3 - 1 ) - %/T53(r 3 ' 3 )l + 0(Pd 3 ) , (89) 



V? a 2 V 21 1 + g r 3 yf] a 2 



15 3(F a 3 ' 3 ) + 0(P(P), 



6tt Pd 
5 r 

r3,l 



2 3(r a 2 fc 2 ) + Q(Z a 2 b 2 ) 



1 pK 1-qPd 2 1 I / 14 st >^2,l^ 
--2 ^ 6V14 9(F afc ' ) - 8d — n{S a ' b ) 



+ 6 3(y^)+3(Z a Y) 



2 V 21 1 + g r 2 a 
%(Y 3 b 3 )+Vi(Z 3 a > 3 )}}+0(Pd 3 ) 



(90) 



(91) 



Eqs. (|8H82p and (f89|) -(|91 |) give the complete harmonic vanishing initial metric perturbations, namely 
decomposition of the initial data (7^ •, if,- ■). We must 



now extract the initial values of the metric perturbations 
and their time derivative by comparing these expressions 
with Eqs. mi)-© and ([7P 1) -([72" 1) . To simplify notation, 
we now drop the truncation error in all equations, since it 
has already been given in the main expansions. Compar- 
ison of Eqs. (|81I82|) with Eqs. (J64J) - (JHBl) yields the non- 



K 2,o 

K 2 
K 3 

K 3,±3 



-2, ±2 



-3,±1 



fK'r = -2 

fh 
fh 



5 (1 + g) 2 ' 
6n g Md 2 1 



1 rr 



5(1 + g) 2 r 3 cr 5 ' 
3tt q(l - g) Md 3 1 



7 3,±3 
rr 



= ± 



7 (1 + g) 3 r 4 a 7 
5nq(l - g) Md 3 1 
7~(l + g) 3 r 4 ~f 



(92) 
(93) 
(94) 
(95) 



As we can see, all the axial metric perturbations ini- 
tially vanish. Now, in order to obtain the initial values of 
the time derivative of the metric perturbations, we must 
compare Eqs. f59" |) -([9"T ]) with Eqs. f70" l) -([72" |) . It is impor- 
tant to realize that in Eqs. (|70 p - (|T2"|) there are terms that 
are associated with the gauge freedom of choosing the 
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slicing, more specifically, terms associated with compo- 
nents of the shift vector [see Eq. Q68p]. Moreover, there 
is no unique way of assigning values to the different time 
derivatives of the metric perturbations and the metric 
perturbations themselves. This reflects the fact that the 
values of the metric perturbations are gauge dependent, 
since in general the components of the metric perturba- 
tions (and their time derivatives) are not gauge invariant. 
Keeping this in mind, we have assigned the following ini- 
tial values to the time derivatives of the metric pertur- 
bations: the non-vanishing polar modes are 



i.2,±2 



■3,±1 



',3, ±3 



2 i 



F Pd 1 

30 ^77 

3tt 1 - q Pd 2 



1 



7 I + q r 4 y/Ja 2 
r»;v l-qPd 2 1 



G 



2, ±2 



K 



2, ±2 



G 



3,±1 




qPd 2 V7 



K 



3,±1 



£3,±3 

and the axial ones are 




1 1 


+ q r 4 a 2 


" 1 - 


-qPd 2 V7 


1 - 


Vq r 4 a 2 


1 - 


qPd 2 ^f 


1 + 


A 9 ' 

q a 


1 - 


qPcPVf . 



q 



= 4V37 — 



•2,±1 



H 



2,±1 




(96) 
(97) 
(98) 
(99) 
(100) 
(101) 
(102) 
(103) 
(104) 
(105) 
(106) 

(107) 
(108) 
(109) 



By using the correspondence between our parameteri- 
zation of the metric perturbations and that of Regge 
Wheeler [Eqs. (|B1[) - (|B5[) in Appendix [B] we have checked 
that up to 0{d 2 ) and O(Pd) our expressions agree with 
those found in [i?! (up to a typo in their hi' ). We have 
decided to assign values to the time derivatives of the 
metric perturbations and the metric perturbations them- 
selves by the following usual convention: all modes with 
£ = 1 are assigned to metric perturbations associated 



with the shift vector. These perturbations represent ei- 
ther translations or rotations of the observers associated 
with the normal to the initial slice with respect to our 
coordinate system. In our case, these perturbations are 
given through the following relationships 



fhtr 



r 



i.±i 



1 + q r 4 a 2 v ' 



Pi 



i.±i' 2 j ±1 ,i,±i . r^-Pd 2 ^ — q 1 

-Pt +h t r =iV67r^ r — -, (111) 

' 1 + q cr 



i.±i 

Pt 



rfh tr 



. f^l-qPd 2 ^J 

' l \hri~, 2 — 2- 112 ) 

V 2 1 + q r A a z 



Here, two comments are in order. First, one can see 
that these equations are consistent in the sense that the 
derivative of Eq. (|112p with respect to r can be reduced to 
a trivial identity by using Eqs. (| 1 1 0[) and (jllip . Second, 
from these equations we can immediately see that the 
shift vector is different from zero [see Eq. {55])]. A non- 
zero shift could in principle be a problem if we wanted to 
place observers at constant r (in the wave zone) , evaluate 
the linear momentum flux, and then infer a recoil velocity 
of the final black hole after the merger. If we were to do 
this, the measured velocity would have a component due 
to the motion of the observers with respect to the position 
of the final black hole, as described by the shift vector. 
This contribution would then have to be subtracted, but 
it can be seen that the shift vector decays quite fast as r 
becomes large and, hence, this effect would be negligible. 



D. Initial data for the master functions 

Using the initial data for the metric perturbations 
[Eqs. (H2)-(IH5l), ([TOI- fTUBl) . and Eqs. fTOTI - fTTO]) ] in the 
master functions [Eqs. (|46p - (|48|) ]. we can compute ini- 
tial data for them: , ¥z) , «™ , , and 

(f£™ M , <j£™ M ) on the initial slice t = t . 

The results for the Regge- Wheeler master function are: 

(113) 



RW " * 



. W£ i^i^^i^Z) (7<r _ 3 ) (U4) 



15 1 + q r 4 



(7 



and for the Cunningham-Price-Moncrief master function 
are 



2,±1 



4 V307T 1-qPd 2 1 



3 5 1 + q r 2 a 



2 ' 



I III 



(115) 
(116) 
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In the same way, the non-vanishing initial data for the 
Zcrilli-Moncrief master functions is given by 



^2,0 

ZM 



Md 2 1 + 5cr 
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3V 5(l + 9) 2 r 2 A 2 cr 5 
Fhv a Md 2 1 + 5cr 
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(117) 
(118) 
(119) 
(120) 

(121) 
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M 



123) 



Note that the master equations do not have the same un- 
controlled remainders as its derivatives, since they come 
from different components of the initial data. In the 
case of unboosted head-on collisions [44[ , the initial data 
scales in powers of d N . Therefore, one only needs to per- 
form one single numerical evolution of the master func- 
tions for some reference value of d — d* , and the results 
for any other value of d can be found using the scaling 
relation. For non-time-symmetric data, such as for qua- 
sicircular or boosted sets, such scaling does not exist. In 
our case, for example, although still scales as d N , its 
time derivative scales as Pd N . Therefore, the mas- 
ter functions themselves are not straightforwardly scal- 
able and several runs with different values of the initial 
parameters must be performed. 



V. RESULTS FROM THE CLA 

In this section, we evolve the master functions with 
the initial data obtained in the previous sections in the 
CLA scheme and report the results for the main phys- 
ical quantities, in particular for the gravitational recoil 
velocities. We first need to choose appropriately the pa- 
rameters that completely determine the initial data, such 
that it describes a binary black hole system merging from 
a quasicircular orbit fsubsection IV A[) . Then, in subsec- 
tion EEI we use a numerical code to evolve the different 
master equations [Eqs. (I4D[) ] and compute the relevant 
physical quantities. We discuss the results and compare 
with previous ones in the literature when possible. 



A. Determining the parameters of the initial data 

Our initial data depends on the following parameters: 

• The total (ADM) mass of the system, M; 

• The mass ratio, where one can use either the bare 
mass ratio q or the physical one Q, related by 
Eqs. dUD and ([15]). 

• The initial separation, where again one can use the 
bare separation d or the physical one D, related by 
Eq. dUD; 

• The linear momentum parameter P of each indi- 
vidual hole. 

Within the family of initial data spanned by these four 
parameters, we need to single out the subset that corre- 
sponds to configurations in quasicircular orbital motion. 
In numerical relativity this is done by looking at the min- 
imum in the binding energy of the system with respect to 
the distance, while keepin g th e total ADM angular mo- 
mentum constant (see, e.g. (52j). We here follow the same 
procedure without using the slow motion approximation. 
The binding energy that we minimize is 



E b = M t 



M 1 -M 2 . 



(124) 



where M. ADM is the total ADM mass and it is computed 
in the asymptotically-flat region containing the two holes 
(S ). This mass is given by (see, e.g. [Zl|l ) 



•Madm = M 



5P 2 

~8/7 ' 



(125) 



where M is given in Eq. (fl"2"|) and /i is the reduced 
bare mass, i.e. /i = m 1 m 2 /M. Moreover, in Eq. (|124|) . 
M 1 and M 2 denote the masses computed in the 
asymptotically-flat regions H 1 and £ 2 . These masses are 
given by (see, e.g. 



Ill) 



Mr 



•Ma 



P 2 

8m. 



(A = 1,2), 



(126) 



where M A is given in Eq. pip . Then, the binding energy 
can be written in the following form 



E h = 



mim 2 



J 



2fid 2 ' 



(127) 



where J is the ADM angular momentum, given by J — 
Pd . This binding energy is formally the same as the one 
corresponding to a binary system in Newtonian gravity. 
One can then minimize this binding energy with respects 
to d, while keeping J fixed, to obtain the condition for 
quasicircular motion (note that in our context there is no 
such a thing as an ISCO) 



J 



H 2 M 



(128) 
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In the same way, one can calculate the associated orbital 
frequency of such orbital motion by differentiating the 
binding energy with respect to J, while keeping d fixed. 
The result is 



J 



(129) 



From Eq. (|128p we can write the linear momentum P is 
terms of the other parameters of our initial data as 



P 



(130) 



The binding energy and other quantities derived from it 
have a Newtonian form because of the particular type 
of initial data that we are using: a conformally-flat 3- 
metric with a Bowen-York extrinsic curvature and a Brill- 
Lindquist conformal factor. The PN metric produced 
by a binar y s ystem differs from conformal flatness at 
0(v A ) (see [72l l73j for the argument in the case of time- 
symmetric initial data), and, hence, the binding energy 
used above differs from the PN binding energy at that 
order. Note, however, that although the binding energy, 
linear and angular momentum used here have a New- 
tonian form, they are not strictly Newtonian. This is 
mainly because the distance parameter d is not the phys- 
ical distance D, which is related to the parameter d via 
Eq. (TTU). 

Adopting Eq. (|130[) for the linear momentum parame- 
ter in our initial data and leaving the total mass M fixed 
(which defines a system of units), we reduce our initial 
parameter space to a 2-dimensional one. The final pa- 
rameter space can be parameterized either by the bare 
quantities (q, d) or by the physical quantities (Q, D) . The 
range of q, or Q, is the obvious one, i.e. [0, 1] , while the 
range for the bare distance parameter d is [d mla , <i CLA ], 
where d CLA is an estimate of the maximum distance for 
which the CLA is expected to be valid. For the case of 
equal-mass head-on collisions, it has been shown [36[, by 
comparing with second-order calculations and with fully 
numerical relativistic simulations, that c! CLA ~ 1.7 M, 
which roughly corresponds to -D CLA ~ 4 M . On the other 
hand, in principle d mln could be just zero however, if we 
adopt the prescription (| 1 30|) for the linear momentum 
parameter, then we are limited by the slow motion ap- 
proximation that we are using, which means that d min 
should be bigger than qM/(l + q) 2 . Finally, we should 
remark that the CLA also is expected to fail in the point- 
particle limit [74j |, but, as we will see, the recoil is very 
small when Q <C 1. 

In order to study the gravitational recoil in the CLA 
scheme, we are going to evaluate the recoil velocity for 
a representative number of (physical) mass ratios for a 
given fixed physical distance D . In particular, we study 
the cases D = 3,3.5,4M and instead of using Q we use 
the physical symmetric mass ratio 



Mi Mo 



V 



Q 



The inverse relation is given by 

Q = i(i-VT^ 



1 . 



(132) 



(Mi + M 2 ) 2 (1 + QY 



(131) 



However, the parameters that appear in our expressions 
for the initial data are the bare ones. Then, in order to 
obtain a plot of the recoil velocity in terms of the physi- 
cal mass ratio, we need to translate from the set (Q,D) 
to (q,d). This, however, is not a trivial calculation be- 
cause the definition of D [Eq. (fl7|) ] is quite intricate, in- 
volving x 1 and x 2 . These numbers are the values of the 
coordinate x in the conformal flat space of the intersec- 
tions of the extremal surfaces (marginally trapped sur- 
faces or apparent horizons depending on the parameters 
of each particular configuration) surrounding each indi- 
vidual hole with the A-axis. The translation has to be 
done numerically through the following iteration scheme 
in which the physical distance D is kept fixed: 

1. Given a value of n we pick an initial guess for the 
bare distance, say d t . 

2. By solving the equations that determine the ex- 
tremal surfaces surrounding each individual holes 
(they are given in Appendix [D]) we find some in- 
tersection points x* and x* 2 . This requires another 
iteration, since we do not know a priori where these 
surfaces are located. What we do is to start, for 
each individual hole, with an initial guess for the 
intersection of the extremal surface at the other 
end of the X axis (the intersection more distance 
to the other hole) and integrate the correspond- 
ing ODEs by using an extrapolation Bulirsch-Stoer 
scheme [7^, [76|, [77}. Then we study whether the 
integration ends far away from the X axis or con- 
verges towards it. We repeat the iteration until we 
find the intersection points x\ and x\ with enough 
accuracy. 

3. Using Eq. (fl~7| we compute the physical distance as- 
sociated with these values of the intersection points 
and (q r ,d^) [where q^ is computed in terms of r\ 
and d. M using expressions (|132[) and fTB])]. . We 
compare and D and stop the iteration if the 
absolute difference between them is smaller than 
10~ 4 M . Otherwise, we go back to point (i) chang- 
ing the ansatz depending on whether is bigger 
or smaller than D . 

We have carried out this iteration for 101 values of 
t] . The coordinate distance (in the conformally related 
fiat space) from the holes to the intersection points of 
the extremal surfaces is shown in Figure [5J Here, we 
observe how these distances move from equal values (top 
right) to the values corresponding to the point particle 
limit, M 2 — > m 2 — » (top left). The bare distance as a 
function of 77 is shown in Figure [3] for the three values of 
fixed proper separation. 
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FIG. 2: Plot of the coordinate distance between the holes and 
the intersection points of the extremal surfaces with the X- 
axis in terms of the symmetric mass ratio r) for three values 
of the physical distance: D = 3 , 3.5 , 4 M . 



FIG. 3: Plot of the bare distance d in terms of the symmetric 
mass ratio r\ for three values of the physical distance: D = 
3,3.5,4M. 



B. Results from the numerical evolution of the 
master equations 

We now have initial data for the master equations and 
also a method to prescribe the initial data parameters in 
a meaningful way. Then, the next step is to evolve the 
master equations [Eq. (|49[) ]. In this paper we use a nu- 
merical code, based on Finite Element methods, that was 
developed in [78| for calculations of the gravitational radi- 
ation emitted by a point particle orbiting a non-rotating 
black hole. This method is based on linear elements and 
hence it has a second order convergence rate with respect 
to the spatial resolution. The time-evolution algorithms 
that it uses are second-order and unconditionally stable, 
since they are based on implicit methods. Apart from the 
tests of the numerical code carried out in [78|, we have 
also done some checks to validate the additional infras- 
tructure added for the gravitational recoil calculations in 
the CLA scheme. First, we have checked that the energy 
and angular momentum emitted in an equal-mass graz- 
ing collision coincide with the ones found by Khanna et 
al in [47j] - Second, we have checked that the recoil veloc- 
ities that we obtain are consistent with the plots shown 
by Andrade and Price [3] for the case of head-on col- 
lisions from rest of unequal mass black holes using BL 
initial data. 

We have then performed evolutions for 101 equally- 
spaced values of the symmetric mass ratio r\ covering the 
whole range [0, 0.25] for the three values of the physical 
distance mentioned above, i.e. D = 3, 3.5, AM . The 



procedure to calculate the bare distance d has been de- 
scribed in the previous subsection. Finally, the linear 
momentum parameter P is obtained through Eq. (|130|) . 
For each evolution we have computed the fluxes of en- 
ergy, angular momentum and linear momentum carried 
by the gravitational waves to infinity. 

Our initial data only has a few non-zero multipoles con- 
tributing to the gravitational radiation emitted. Thus, 
the expressions for the different fluxes simplify dramati- 
cally. The energy flux is given by 
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8^ 
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(133) 



Figure [5] shows the total energy to infinity, given by the 
integral of Eq. (|133p over time, as a function of the sym- 
metric mass ratio. 

The angular momentum flux also simplifies greatly and 
becomes 



5 
2 

15 r 

~2 



Note that Eq. (|134[) does not contain any contributions 
from the axial modes, since the only non-zero axial mode, 
$cpm, is purely real. Figure [5] shows the total angular 
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FIG. 4: Energy radiated to infinity in terms of the symmetric FIG. 5: Angular momentum radiated to infinity in terms of 



mass ratio r\ for three values of the physical distance: D 
3,3.5,4M. 



the symmetric mass ratio rj for three values of the physical 
distance: D = 3 , 3.5 , 4M . 



Eqs. (|57H59|) . These expressions reduce to 



momentum radiated to infinity, given by the integral of 
Eq. (|134p . as a function of the symmetric mass ratio. 

Finally, the gravitational waveform also simplifies, and 
we obtain 



+ 2R(*|=) 5i(- 2 y 2 ' 2 ) - 23(4< 2 i) 3(_ 2 y 2 > 2 )] 

+ 5i(* 3 M 3 ) R(- 2 K 3 ' 3 ) - 3(- 2 r 3 ^ 3 )] } (135) 
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where the definition and some properties of the spin- 
weighted spherical harmonics g Y m are given in Ap- 
pendix IA 31 Note that the x -polarization consists purely 
of the axial modes, while the +-polarization contains 
only polar contributions. Figure [B] shows a typical met- 
ric waveform, namely h + as a function of time, for an 
observer located at ~ 300M on the Z-axis (the cross po- 
larization vanishes on this axis). 

Let us now concentrate on the main physical quantity 
of interest, namely the linear momentum flux, given in 
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(139) 



As we can see, there is only contributions from the over- 
lap of polar modes with different t and m. From this 
flux, the recoil velocity can be obtained by performing 
the following integration 



(140) 
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FIG. 6: Metric waveform h + for the case D = 3.5M and 
n — 0.185 as a function of time. The observer is located at 
~ 300M on the Z-axis 



where the integration times, t i and tp are such that the 
time interval includes essentially all the contribution from 
the waves to the flux. We can then calculate the magni- 
tude of the recoil velocity simply by 

v xecoil = ^(v x ) 2 + <y) 2 + (v z ) 2 , (141) 

where v z = in our case, due to the choices made in the 
initial setup. Figure [7] shows the time derivatives of the 
master functions that contribute to the recoil velocity for 
a typical evolution. In this figure, we have separated the 
real (bottom panel) and the imaginary (top panel) parts 
of these time derivatives. Observe that the magnitude 
of the £ = 2 modes is much bigger than the one of the 
£ = 3 modes, as expected. This also gives an indication 
that the superposition of the I = 2 modes with £ — 2 
and £ = 3 modes is going to be the dominant contribu- 
tion to the gravitational recoil. The contribution from 
superpositions involving higher I 's is going to be much 
smaller. 

Figure [8] shows the the linear momentum flux as a 
function of time. Observe that the magnitude of the 
x-component is bigger than the ^-component, which re- 
flects the fact that our configuration corresponds to the 
transition from merger to plunge. 

Finally, Figure O presents the magnitude of the recoil 
velocity as a function of the symmetric mass ratio, for the 
following initial physical separations: D = 3 , 3.5 , 4M. 
For all cases studied, the maximum velocity is reached for 
a symmetric mass ratio of n ~ 0.19, which agrees with 
the value reported in Refs. [12, [H| up to uncontrolled 
remainders. Observe that this maximum is not a strong 



FIG. 7: Time derivate of the Zerilli-Moncrief and 
Cunningham-Price-Moncrief master functions as a function 
of time, for the case D — 3.5M and r\ = 0.185. The plots in 
the top panel represent the imaginary parts whereas the ones 
in the bottom panel represent the real parts. 



peak, but instead resembles a plateau, where this maxi- 
mum is spread out for a number of etas, as also seen in 
other calculations [H, . 

VI. ESTIMATING THE TOTAL RECOIL 

In this section, we discuss the recoil velocities obtained 
from the evolution of the master functions and produce 
lower and upper limits for the total recoil velocity. In 
particular, we will provide analytic approximations to 
the data and we will also compare these results to other 
ones already present in the literature. 

One of the limitations, and at the same time an ad- 
vantage, of the CLA scheme is that the initial separa- 
tion of the black holes must be sufficiently small in some 
well-defined sense. Apart from numerical relativity, this 
method is the only known one to be capable of produc- 
ing accurate estimates of physical quantities near plunge. 
This advantage, however, is a double-edged sword since 
the method cannot account for the inspiral phase. Ac- 
tually, the initial separation must even be smaller than 
that for which the last ISCO exists. Thus, not only is 
the inspiral phase neglected but also the beginning part 
of the merger phase. 

Due to these limitations, an approximate value for the 
total recoil velocity cannot be provided by the CLA alone, 
without supplementing it with some other scheme valid 
when the system is well separated. The PN scheme is 
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well suited for this task and extensive studies have been 
recently carried out [H, [23| to estimate the recoil veloc- 
ity. The approximate recoil velocity accumulated from 
infinity up to some final separation in the PN scheme is 
given by [22J 



FIG. 8: Linear momentum fluxes, P<? w and Pq W , as a function 
of time, for the case D = 3.5Af and rj = 0.185 . 
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with remainders of 0(v 5 ). In Eq. Ijl42|) . Xf = (muf)" 2 
is a PN parameter, m PN = jt^ pn + m 2 PN is the to- 
tal mass, m 1 2 PN are the masses of the PN point par- 
ticles, r\ = TTi\ui2 1 TTi is the symmetric mass ratio and 
<5to pn = ra 1 PN — m 2 PN is the mass difference. The PN 
masses mi ; 2, PN have been shown to agree, within the PN 
approximation, with the horizon masses M\ i2 [z2,[73| and 
we make this identification here. The angular velocity to 
is given to 0(v 4 ) by 



where v PN [D 2 , oo 



} is the PN estimate for the recoil ve- 
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(143) 

and uif is the angular velocity evaluated at some final co- 
ordinate separation b / . Post-Newtonian theory is usually 
carried out in harmonic coordinates, which are different 
from the Schwarzschild coordinate system we use in the 
CLA scheme. However, sufficiently far from the holes, 
D ~ b, to 0(v 2 ). 

Supplementing the CLA estimate with the PN esti- 
mate, we can obtain upper and lower limits on the possi- 
ble values of the magnitude of the recoil velocity. A lower 
limit can be obtained via 



u low = «c LA [0,4M]+ UpN [6M,oo] 



(144) 



locity of Eq. (IL42[) evaluated at bf = D 2 - For this lower 
limit, we evaluate the PN estimate at the edge of the re- 
gion of validity of the PN approximation, i.e. bf = 6M, 
or equivalently Xt = 6 -3 / 2 , as done in Ref. (22J. This lo- 
cation corresponds to the ISCO of a test particle around 
a Schwarzschild hole of mass M. One obtains this value 
of Xf(bf = 6M) by neglecting terms of 0(v 2 ) and higher 
in Eq. (|L43|) . If we had included these higher order terms 
in luj and Xf, the upper bounds would have decreased 
by approximately 50 km/s. These higher order terms, 
however, become large as b becomes smaller, and thus 
we choose to neglect them to have a conservative upper 
bound. In Eq. (|144|) . -y CLA [0,D 1 ] is the estimate of the 
recoil velocity in the CLA approximation with an initial 
proper separation of D = D 1 . 



The estimate of i> low is a lower limit because it does not 
take into account the contribution to the gravitational 
recoil in the region b € (4, Q)M. In this region neither 
the CLA, nor the PN scheme, is guaranteed to provide an 
accurate estimate for the recoil. However, it is possible 
to construct upper limits by modelling either the entire 
region or part of it with PN and CLA estimates. Such 
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FIG. 10: Estimated lower (black solid curve) and upper lim- 
its 1 (dotted red line) and 2 (dashed blue line) for the recoil 
velocity after a binary black hole merger as a function of the 
physical symmetric mass ratio. Note that the maximum oc- 
curs roughly in the same place, namely rf ~ {0.19, 0.2}. 



FIG. 9: Magnitude of the recoil velocity in terms of the sym- 
metric mass ratio r\ for three values of the physical distance: 
D = 3,3.5,4M. 



upper limits are given by 

V! = «c LA [0,4M]+ v PH [4M,oo], (145) 
K P , 2 = « CLA [0,5M]+« PN [5M,oo]. (146) 

These expressions are upper limits because the contri- 
bution to the recoil estimated either with PN theory or 
the CLA approximation in the region &/ S (4,6)M is 
monotonically increasing with bf. 

Equations (fl"44l) , ([145]) . and (fi"46"l) require some extra 
justification and clarification. In general, it is not true 
that the magnitude of the total recoil can be estimated by 
adding the magnitude of the integrated momentum flux 
in the region [4M, oo] to that in the region [0,4M]. The 
important point is that the cut is made at a separation 
D — AM in the regime where the main contribution to 
the recoil comes from. Then, the main contribution to 
each recoil velocity vector comes from the region near the 
cut and hence, the error we made by adding the norms 
will be relatively small. Independent of this argument we 
have the inequality i>[0, oo] < v[0, D cut ] +v[D cllt , oo] (with 
equality when v[0, D cut ] and v[D cut ,oo] are aligned). In 
this sense, the proposed upper limit is indeed always an 
upper limit, irrespective of the orientation of the vectors. 
As for the lower limit, neglecting the accumulated recoil 
in the region [4M, 6M] is a very conservative estimate, 
because there the recoil accumulates greatly. Thus, the 
issue of the orientation of the vectors will not affect the 
fact that this is a lower limit, as other recent estimates 
in the literature confirm. 

Figure [TU] shows the behavior of these upper limits (red 
dotted and blue dashed lines respectively) and the lower 



TABLE I: Values of the parameter of the non- linear fitting for 
the following models: the CLA with initial separations of D = 
{3, 3.5, 4}M; the lower and upper limits of Eqs. (TTJU) -(TJJiJ; 
Taylor PN (BQW) and EOB PN (DG) calculations [22, 23]. 



Model 


a (km/s) 


b 


c 


Mean square error 


« CLA [0,3M] 


1841 


-3.31 


3.45 


0.001 


« CLA [0,3.5M] 


3548 


-3.15 


3.33 


0.003 


« CLA [0,4M] 


6576 


-2.98 


3.21 


0.008 


*>i<™ 


7782 


-2.51 


2.73 


0.008 


«up,l 


14802 


-1.13 


1.48 


0.008 


«up,2 


23124 


-2.33 


2.61 


0.06 


V BQW 


12891 


0.25 





10~ 8 


VDG 


4483 


-0.95 


2.68 


io- 10 



limit (black solid curve as a function of r\. The maxi- 
mum in these curves occurs roughly at the same sym- 
metric mass ratio, namely r\ ~ {0.19,0.2}. The slight 
disagreement in this maximum is within error bars and 
rooted in that PN theory predicts it at approximately 
r\ — 0.2, while the CLA predicts it at 77 — 0.19. We 
should note that the maximum recoil from u CLA [0,4M] 
and i> CLA [0, 5M] is approximately 64 km/s and 215 km/s, 
while the maximum recoil from v PN [4M, 00] , v PN [5M, 00] 
and f PN [6M, 00] is approximately 160 km/s, 50 km/s, 
and 20 km/s respectively. 

A non-linear fit can be performed to these curves via 
Eq. Q) 

v m = arf \J\ - Ar) (l + brj + cr] 2 ) , 

where the fitting parameters a, b and c are listed in Ta- 
ble D 

Observe that the mean square error for all cases is 
small, which is an indication that Eq. (JTJ is a good ana- 
lytic model for the functional form of the recoil velocity. 
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In this table, we also present the values corresponding to 
the estimates of Refs. [3 (BQW) and [H (DG.) Since 
the predictions of these references are based on analytic 
formulae, the mean square error can be made arbitrarily 
small by increasing the number of points in the discretiza- 
tion of the analytic curve. 

With the analytic fits to the upper and lower limits, 
we can construct a curve that is in between these limits 
with an error given by the distance from the curve to the 
upper or lower bound. Such a curve is given by Eq. ([1]) 
with the following fitting parameters 



300 r 



a u „b u 



a,„ w c, f 



(147) 
(148) 
(149) 



while the error on this curve is also given by Eq. fl} with 
the following fitting parameters 



(150) 
(151) 
(152) 



a,„ w c,. 



■ a„„c„ 



where the subscript low and up stand for the fitting pa- 
rameters of the lower or upper limit respectively. This 
curve is only an alternative way to visualize the upper 
and lower limits of Fig. [TO] The curve is not to be inter- 
preted as the best guess in this work, since in principle, 
the recoil velocities present in nature could be closer to 
either upper or lower limit. 

We can now compare these estimates of the recoil ve- 
locity with those present in the literature. Fig. QT] shows 
the recoil velocity in units of km/s as a function of the 
physical symmetric mass ratio as estimated in the litera- 
ture and by this paper. As is clear from the figure, there 
are many approaches to calculate this velocity and not 
all of them agree. The symbols used are the following: 
squares and triangles stand for the results obtained us- 
ing black hole perturbation theory in the extreme-mass 
ratio approximation in Refs. [TtI [l3| respectively; circles 
stand for the calculations carried out via the Lazarus ap- 
proach [25[> stars and crosses correspond to the results 
coming from a full numerical relativistic simulation (PSU 
stands for Ref. [Hj| and NASA stands for Ref. H); the 
dotted line and the dashed line correspond to the 2 PN 
Taylor expansion approach [22] and the 2 PN effective- 
one-body (EOB) approach [23j respectively. The solid 
line with error bars is the estimate of Eq. (|149I152[) that 
properly condenses the lower and upper limits into one 
curve. We briefly describe each a ppr oach below. 

In the PN calculations of Ref. (22j , the recoil velocity 
(dotted line in the figure) and momentum flux are esti- 
mated by studying the 2 PN Taylor-expanded radiative 
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FIG. 11: Estimates for the recoil velocity (km/s) of the in- 
spiral and merger of a binary system of compact object as 
a function of the physical symmetric mass ratio parameter. 
The symbols used are the following: squares and triangles 
stand for black hole perturbation theory results [l7l , Il9j ; cir- 
cles stand for Lazarus results [251 ] : stars and crosses corre- 
spond to full numerical relativistic simulation [3l], H^; the 
dotted line and the dashed line correspond to 2 PN Taylor 
expansions [13] and 2 PN effective-one-body expansions [23|] 
respectively. The solid black line corresponds to the estimate 
of this paper, which, together with the error bars, condense 
both upper and lower limits. Other error bars, when present, 
correspond to an estimate of some of the error contained in 
the calculation. 



moments of a binary system of compact objects, while in 
Ref. [HI an effective-one-body approach is used (dashed 
line in the figure.) Post-Newtonian calculations are usu- 
ally valid only when the binary is weakly gravitating, or 
equivalently when the orbital separation is greater than 
the ISCO. In this regime, the recoil velocity (Eq. (|142[) 
has been found to be small for any mass ratio (usually less 
than 20 km/s), since, as expected, most of the contribu- 
tion to the recoil comes from the merger part of the inspi- 
ral. In Ref. [13] the calculation is extended through the 
merger by integrating the 2 PN Taylor-expanded momen- 
tum flux along a geodesic of the Schwarzschild metric. 
On the other hand, Ref. [23j uses the effective-one-body 
Hamiltonian to extend the inspiral through the merger. 
Both of these approaches have inherent errors that are 
difficult to estimate without calculating the 3 PN contri- 
butions to the recoil velocity. 

Black hole perturbation theory has also been used to 
estimate the recoil velocity in Refs. [13, EH- In these 
studies, the extreme-mass ratio approximation is adopted 
(i.e., Q <C 1) and then the system is approximated as a 
point particle orbiting a black hole. The first study of 
the recoil velocity using this formalism was performed 
in Ref. [lj| (squares in the figure), but there the grav- 
itational force was treated as Newtonian _and only the 
lowest multipoles were considered. 



In Ref. [17| , these rel- 
ativistic effects were taken into account, as well as spin, 
and the velocity estimates were improved (triangles in the 
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figure.) The extreme-mass ratio approximation, however, 
requires Q <C 1, which allows the exploration of a limited 
section of the ?7-space. 

A combination of black hole perturbation theory and 
full numerical relativity (the so-called Lazarus approach) 
has also been implemented to estimate the recoil veloc- 
ity [2||. In this case, a full numerical relativistic simu- 
lation is carried out until the black holes merge and a 
single perturbed spinning black hole has formed. Then, 
this spacetime is used as initial data in a Teukolsky evo- 
lution to determine the recoil velocity (circles in the fig- 
ure.) The error in this calculation is rooted in the inter- 
pretation of the initial data as that of a perturbed Kerr 
spacetime. Finally, there have also been recently full nu- 
merical relativistic simulations of binary black hole coa- 
lescence jU, H2| (shown as stars and crosses respectively 
in the figure.) In this case, the error shown in the figure 
is assumed to be given only by finite differencing, while 
the error due to initial data is neglected. 

Even though there has been much work in the calcula- 
tion of the recoil velocity there is still some disagreement. 
In Fig. [TT] we observe that there seem to be three groups 
of results: one that clusters around the 2 PN Taylor ex- 
panded result; another that is close to the 2 PN effective- 
one-body result; and a third one that is in between the 
first two. This disagreement, however, is misleading in 
several ways. First, some estimates of the recoil velocity 
quote no error bars, as is the case of the first pertur- 
bation theory approach [l9j and the 2 PN effective-one- 
body approach [231 ] . Second, the error bars that do exist 
in other calculation are only estimates and could very 
well have been underestimated. A surprising disagree- 
ment is between the PN approaches, since when used to 
calculate other quantities, such as the angular frequency 
at the ISCO, they do agree. This disagreement seems 
to be rooted in the fact that the greatest contribution 
to the recoil velocity comes from the merger part of the 
inspiral, where neither extension of the PN approach is 
guaranteed to be accurate. 

Our estimates seem to agree with most results if one 
accounts for error bars. The results of Refs. 
(squares, triangles and circles in the figure) seem to over- 
estimate the recoil, which is expected in the case of the 
extreme-mass ratio approximation. The PN results of 
Refs. [iH [23| seem to overestimate and underestimate 
the recoil respectively, but they are consistent with our 
bounds if one takes their error bars (not shown in the 
figure) into account. The full numerical relativity results 
seem to overlap with our bounds, although there arc only 
a few of them. We should note that our bound seems 
to disagree with the full numerical relativistic result for 
r\ ~ 0.18, but that result seems to be an underestimate 
because of the small initial separation [79l |. 



VII. CONCLUSIONS AND DISCUSSION 

We have calculated the recoil velocity after the merger 
of an unequal mass binary black hole system using the 
CLA scheme. This approximation assumes that the black 
holes are close enough that the system can be approxi- 
mated by a single perturbed black hole spacetime. In 
contrast to other approaches, except for full numerical 
relativity, this approximation allows us to make valid 
statement about physical process when the system is 
close to plunge. Therefore, it is of great interest to use 
this method for the study of gravitational recoil. How- 
ever, the CLA has the disadvantage that it cannot be 
used during the beginning of the merger or the inspiral 
phases. 

Initial data for the CLA can be constructed analyti- 
cally by mapping data suitable for a binary black hole in- 
spiral to that of a single perturbed hole. With such initial 
data, the Cunningham-Price-Moncrief and the Zerilli- 
Moncrief master functions can be numerically evolved 
from some initial proper separation through ringdown. 
These gauge-invariant master functions contain all the 
information necessary to evaluate the gravitational met- 
ric waveforms and, thus, the energy, angular momentum 
and linear momentum fluxes carried away from the sys- 
tem. 

The results obtained can be summarized as follows. 
First, the maximum recoil velocity obtained in the CLA 
scheme is of v ~ 64 km/s for the maximum initial sep- 
aration allowed by this method (D = AM). This max- 
imum occurs at a symmetric mass ratio of 77 ~ 0.19. 
By supplementing this estimate with PN ones valid in 
the inspiral regime, we obtain lower and upper bounds 
with maxima of v ~ 84 km/s and v ~ 220 km/s respec- 
tively. We have further provided non-linear analytic fit 
functions that conveniently parameterize these bounds, 
together with the results from the CLA, and other re- 
sults in the literature. These results also suggest that 
there is a region around D ~ {4, 6} M that greatly con- 
tributes to the recoil, but can only be poorly modelled 
by current approximation schemes. 

Ultimately, the estimates presented here suffer of the 
same predicaments as other calculations. Due to its ana- 
lytical nature, the CLA relies on certain assumptions that 
do not hold over the entire history of the binary. Such 
assumptions introduce an error in the estimated recoil 
that is difficult to quantify. In particular, the assump- 
tions made here are the following: close separations ; 
slow-motion; simple initial data. The close-limit assump- 
tion is essential to allow a mapping of a binary inspiral 
to a single perturbed spacetime. The slow-motion ap- 
proximation supplements the close-limit assumption and 
can, in principle, be improved on in future extensions 
of this work. The choice of initial data is assumed to 
represent the gravitational content of some initial slice, 
although we know that this fails even at large separation 
because it does not agree with the deviations from confor- 
mal flatness predicted by PN theory. Moreover, it does 
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not contain any radiation, which is not what it should 
be expected for initial data corresponding to a snapshot 
of the orbital evolution. Due to these assumptions, the 
estimate of the recoil velocity will be contaminated by 
some error. However, experience in CLA calculations in- 
dicates that the error made only overestimates the phys- 
ical quantities calculated, relative to full numerical sim- 
ulations m, S3, 113. 

Future work will concentrate on extending this ap- 
proach to second order in P and to other, more real- 
istic, initial data sets [U [Hj]. Ultimately, it would be 
interesting to compare the CLA approach directly to full 
numerical relativistic simulations in an attempt to deter- 
mine the region of validity of the CLA more accurately. 
Another possibility is to use a multi-parameter pertur- 
bation scheme (see [H, H3, |H, HH) where perturbations 
in the linear momentum and separation parameters can 
be cleanly separated at the different perturbative orders. 
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1. Special polynomials 

The expression for the associated Legendre polynomi- 
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where I is a non-negative integer and m is an inte- 
ger restricted to the following range: m S (—£,—(£ — 

,e-i,e). 

Gegenbauer polynomials, also known as ultraspherical 
harmonics [57j, are generalization of Legendre polyno- 
mial for higher dimensional spaces. They can be written 
in terms of other special functions, as in 

w = F(A + l/2) F(n + 2A) (A _ 1/2 , A _ 1/2) 
r(2A) F(n + A + l/2) n 

where A is a real number, n is a positive integer, Y is 
the Gamma function, and P^ 1 ^ 2 ' are the Jacobi poly- 
nomials • There are also recursion relations for these 
polynomials, but we will not present them here. Instead, 
we will provide the first few Gegenbauer polynomials 
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(A4) 
(A5) 
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These polynomials allow for the far field expansion of 
potentials that scale as \x — xo\~ a to arbitrary order. 



2. Spherical Harmonics 

The scalar spherical harmonics are solutions of the 
eigenvalue problem described by the following equation 



Note added at revision: After this paper was sub- 
mitted, full numerical relativity estimates have ap- 
peared 88] , and a best-bet estimate based on our scheme, 
which includes some refinements in the PN part with re- 
spect to the present paper, has been compared to them 
in [8^ ]. The agreement between them is remarkable and 
supports the approach presented in this paper. 
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(A7) 



where (£, m) have the same range of values as in the asso- 
ciated Legendre polynomials above. In this paper we use 
the conventions of [77], l9fj to define specify the solutions 
of Eq. (|A7[) . The precise expression is given by 
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APPENDIX A: CONVENTIONS FOR SPECIAL 
FUNCTIONS 



In this appendix, we describe the conventions we use 
for the special functions presented in this paper, and we 
also present some important properties of such functions. 



The scalar spherical harmonics form an orthonormal ba- 
sis on the two-sphere, that is 



(A9) 



where 5 ab denotes the Kronecker delta. The vector spher- 
ical harmonics are defined in terms of the scalar ones as 
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in Eq. (|37p . and from this definition we can derive the 
following orthogonality relations: 



dnn ab Y* m Yf m ' = £(£ + l)S u 'S mm ' 
s 2 

dtiSl ab Si m Si' m ' = £(£ + 1)6"' 6 mm ' 



dnn ab Y^ m s e b ' r 

s 2 



= 0. 



(A10) 
(AH) 
(A12) 



The (symmetric) tensor spherical harmonics used in this 
paper are also constructed from the scalar ones by means 
of Eqs. (f38|) and (f39|) . from where the following orthogo- 
nality relations can be deduced: 
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and 
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3. Spin-weighted scalar spherical harmonics 

Spin-weighted scalar spherical harmonics are another 
basis to expand functions on the 2-sphere. They can be 
defined by the following general formula 91] 
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It is important to mention that these definitions are ap- 
plicable only to integer powers of the spin- weight. A 
generalized definition for half-integer powers of s exists 
but will not be discussed here (see [911]). 

In this paper we are interested in the s = — 2 case, 
for which the definition of the spin-weighted spherical 
harmonics reduces to [93l 
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APPENDIX B: RELATIONS WITH THE 
REGGE- WHEELER PARAMETERIZATION OF 
THE PERTURBATIONS 

For the sake of completeness, we give here the rela- 
tions between our parameterization of the metric pertur- 
bations and the one used by Regge and Wheeler [58| • For 
polar modes (our notation is on the left column and the 
one of Regge and Wheeler is on the right one) we have 
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where d (d) is a ladder operator, usually called the edth 
operator, that raises (lowers) in a unity the spin weight 
of any quantity. Its action on a scalar Q can be expressed 
in the following way 



dQ = m a d a Q- -(m a m b V b m a )Q 1 (A19) 
dQ = m a d a Q + ^{m a m b V b rh a )Q , (A20) 

where s is the spin weight of Q and (m a , rh a ) is a 
null complex basis on the 2-sphere {VL ab m a m b = , 
£l ab m a "fh b = 1). It is worth noting that the action of the 
edth depends explicitly on the spin weight of the quan- 
tity on which it acts. Taking m a — i [l, -^-g] , we can 
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The expressions for the master functions, in 
Schwarzschild coordinates, in terms of the parameteri- 
zation of Regge and Wheeler are: 
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APPENDIX C: ON THE DERIVATION OF THE 
LINEAR MOMENTUM FLUX FORMULA 

In order to obtain Eqs. ([57]) -([59 ]) from Eq. ([55]) we need 
to use the decompositions of products of spherical har- 
monics in single harmonics typical of problems that deal 
with angular momentum coupling (for accounts dealing 
with this problem see [67| , where multipole expansions of 

I 



At this point everything reduces to evaluating the inte- 
grals on the 2-sphere. To that end, the starting point is 

I 



where the objects with the round brackets are the 
3j-Wigner symbols, which are related to the Clebsch- 
Gordon coefficients. They are subject to certain selec- 
tion rules, namely, (£,m) , (£',m'), and (L,M) are inte- 

I 



where C(£,£',L) is a constant given by 

C(£, £', L) = ^{L 2 {L + l) 2 +£ 2 (£ + l) 2 + £' 2 {£' + l) 2 
+ 2L(L + 1) - 2 [£(£ + 1) + £'(£' + 1)] [L(L + 1) + 1]XC4) 
A similar formula can be found for the second term 



gravitational radiation in different sets of harmonics are 
described; see j94| for a recent systematic treatment of 
higher-order perturbation theory where these issues are 
also treated). 

Introducing Eq. (f5"l)) into Eq. (|5"5")1 , using that any 
spherical harmonic S lrn has the property S e, ~ m = 
(— l) m S m , and the fact that Z^™ is symmetric and 
traceless we get 



(CI) 



I 

the well-known formula 19 



(C2) 



I 

gers with the usual ranges of values; m + ml + M = ; 
and the triangular inequality \£ — £'\ < L < £ + £' . By 
using ([C2j and the definition of Z*f [Eq. (38])] we find 
the following relationship 
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I 

in (|C1|) . Then, since the components of r^ bs are linear 
in Y 1 '" 1 , the integral that appears in Eq. (|C1[) is now 
straightforward. In order to get Eqs. (|57 |) -([59 | we just 
need to use the selection rules of the 3 — j Wigner sym- 
bols and the following additional properties: 
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where J = ji + j 2 + j • 

Alternatively, one can use the following relation for the 
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integral of three spin- weighted spherical harmonics: 
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APPENDIX D: DETERMINING EXTREMAL 
SURFACES IN THE BRILL-LINDQUIST BINARY 
BLACK HOLE DATA 

We are interested in finding the location of the ex- 
tremal surfaces that surround the individual holes in the 
BL binary black hole initial data that we are using in this 
paper. When the holes are separated enough these sur- 
faces form the apparent horizon of the initial data (since 
we are neglecting the extrinsic curvature the data is time 
symmetric). If we put the two holes close enough, an- 
other maximal surface enclosing the two holes appears 
and becomes the apparent horizon, and then, the two in- 
dividual maximal surfaces are called marginally trapped 
surfaces. 

Following Bishop (9(| [97| , in order to look for maximal 
surfaces in the BL data it is very convenient to exploit the 
cylindrical symmetry of the configuration by expressing 
the metric in cylindrical coordinates (for coherence with 
the conventions of the paper we choose the axis of sym- 
metry to be the X axis) : x — x , y = p cos $ , z — p sin i? . 
Then, the line element of Eq. @ becomes: 



ds 2 
where now 
$ = 1 



^{p.x) {dp 2 + p 2 dd 2 + dx 2 ) 



2v / p 2 + {x _ xj2 2 Vp 2 + {x - x 2 y 



(Dl) 
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reduces to that of finding a path (p(A) , x(X)) in the sub- 
space (p , x). The area of a surface with cylindrical sym- 
metry can be written as 



A= 2tt 



d\, 



(D3) 



where the dots denote differentiation along the path, that 
is d/dX, and (X 1 , X 2 ) are the intersections of the surface 
with the symmetry axis X . The equations for the path 
(p(A) , x(XY) is found by extremizing the area, 5A = , 
Bishop [9a, took A to be an affine parameter, that 
is, such that A = 2n (A 2 — A x ) . We have fixed A in such 
a way the ordinary differential equations for (p(A) , a;(A)) 
are simple and amenable for numerical computations. In 
particular, we have chosen A so that 



V ' P 2 + x 2 = p $ 4 , 
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which is a constraint preserved by the Euler-Lagrange 
equations that we obtain from 5 A = 0: 



-d p V, 



where the potential V is given by 

V(p,x) = ±p 2 &. 
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Moreover, the problem of finding the extremal surfaces 
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